function rho = inprop(om,gam,nd,t,rho,obasis,coef)
% Propagation of density matrix with dissipation! 
%  currently onlu single step
% The program will be as general as possible, but in will be 
% devoted to calculate the simplest case in Huisinga article. 

nleja=10;
ntotal=10*nleja;
% nleja number of interpolation (Leja) points.
% ntotal number of points to generate Leja points.
% zl the Leja points.
% zdvd the divided differences.

      r=2.d0*om/gam;

      r=(1.d0-r)/(1.d0+r);    
% Reescaling the Liouvillian.
% L -----> L/sc, shuch as the spectrum is into the scaled domain.
% In this case is very simple because we know how is the spectrum (notes).


%Warning! for space representation it will be necessary to stimate the
% most energitical state represented.

      sc=((nd-1.)*om)/(1.d0-r);

% Reescaling the Liouvillian... that is equiv to rescale the coef.

      omsc=om/sc ;
      gamsc=gam/sc;

% Reescaling time step      
      tsc=t*sc;

% Generating Leja points (the algorithm depends on the 
% elliptical shape choosen; subroutines cpoints,zshape).       
      [zl,lc]= leel(ntotal,nleja,r);
      zdvd = dvd(nleja,zl,tsc);

%  Time propagation. (time step is implicit in zdvd). CCCCCC
   rho = dpropa(nd,rho,omsc,gamsc,nleja,zl,zdvd,obasis,coef);
% Output: rho(t1+t)
      end 
